subroutine faraday
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
	real(double) :: dummy
!
!     ******************************************************
!
!     solve Faraday's law to advance the magnetic field (cell centered)
!
         call curlc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
            , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, vol, ex, ey, ez, curlx, curly, curlz)
!
!
!	apply bc on B cell centered
!
!	Note that the last entry is set to 0 to impose zero tangential value
!       in the ghost and boundary centres, 1 to apply neumann there. Zero value needs
!	to be applied to the curl E but 1d0 needs to be applied to full B
!


      call boundary_b_centered(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         curlx, curly, curlz, periodicx,0d0)
!
         dtc = clite*dt
!
!         do n = 1, ncells
!
!            ijk = ijkcell(n)
!
         do i=1,ibp2
         do j=1,jbp2
         do k=1,kbp2

            ijk= 1+ (i-1)*iwid+(j-1)*jwid+(k-1)*kwid

            bxn(ijk) = bxn(ijk) - curlx(ijk)*dtc
            byn(ijk) = byn(ijk) - curly(ijk)*dtc
            bzn(ijk) = bzn(ijk) - curlz(ijk)*dtc
!
         end do
         end do
         end do
!
!     ****************************************************
!
!	apply bc on B cell centered
!
      call boundary_b_centered(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxn, byn, bzn, periodicx,1d0)
!
      call vtxb (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, vol, bxn, byn, &
         bzn, bxv, byv, bzv)
!
      call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, wate&
         (1,22))
!
      call bcperv (ibp2, jbp2, kbp2, iwid, jwid, kwid, sdlt, cdlt, bxv, byv, &
         bzv, periodicx)

      call boundary_b_vertex(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxv, byv, bzv, periodicx,1d0)
!
      call current

!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         jx0, jy0, jz0, periodicx)
!
!
      call bdrift (ncells, ijkcell, nsp, itdim, iwid, jwid, kwid, qdnv, qom, &
         bxn, byn, bzn, dt, clite, wate(1,1), wate(1,2), wate(1,3))
!
!     impose toroidal boundary conditions on "alfab"
!
      dummy = 0.0
!
      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, wate(1,1), wate(1,&
         2), wate(1,3), periodicx)
!
      call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, wate(1,1), wate(1,2), wate(1,3), udrft, vdrft, &
         wdrft)
!
      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, udrft, vdrft, &
         wdrft, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, udrft, vdrft, wdrft)
!
      write (*, *) 'vinit:renormalize udrift'
      do n = 1, nvtx
         ijk = ijkvtx(n)
!
         rvvol = 1./vvol(ijk)
!
         udrft(ijk) = udrft(ijk)*rvvol
         vdrft(ijk) = vdrft(ijk)*rvvol
         wdrft(ijk) = wdrft(ijk)*rvvol
!
      end do
!
      return
      end subroutine faraday


subroutine div_cleaning
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
	real(double) :: rcntr, turnon
!
!
!     advance electric field from E(t+cntr*dt) to E(t+dt)
!
      rcntr = 1./cntr
      if (eandm) then
         turnon = 1.
      else
         turnon = 0.
      endif
!
!
      do n = 1, nvtx
         ijk = ijkvtx(n)
!
         ex(ijk) = (ex0(ijk)+(ex(ijk)-ex0(ijk))*rcntr)*turnon + (1. - turnon)*&
            ex(ijk)
         ey(ijk) = (ey0(ijk)+(ey(ijk)-ey0(ijk))*rcntr)*turnon + (1. - turnon)*&
            ey(ijk)
         ez(ijk) = (ez0(ijk)+(ez(ijk)-ez0(ijk))*rcntr)*turnon + (1. - turnon)*&
            ez(ijk)
!
         ex0(ijk) = ex(ijk)
         ey0(ijk) = ey(ijk)
         ez0(ijk) = ez(ijk)
!
      end do

!     calculate the source for Poisson's equation
!     the call computes a correction to E0
!
               call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
	dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk) - divpix(ijk)
		  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               enddo
	       write(*,*)'prima=',dnorm
!
!     with these parameters, poisson solves poisson's equation
!
               dumdt = 0.0
               dumscal = 1.
               wk = 0.
               phibc = 0.
!
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,1), wate(1,2)&
                  , wate(1,7), wate(1,8), itmax, error, wate(1,10), dumdt, &
                  dumscal, wk, phibc, ex, ey, ez, wate(1,4), wate(1,5), wate(1,&
                  6), phipot, wate(1,9), wkl, wkr, wkf, wke, periodicx)
!
!
               call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, &
                  c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
                  c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, phipot, gradxb, &
                  gradyb, gradzb)
!
!
               call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxb, gradyb, &
                  gradzb) !....nothing is done with this call
!
!dir$ ivdep
               do n = 1, nvtx
                  ijk = ijkvtx(n)
                  rvvol = 1./vvol(ijk)
                  ex0(ijk) = ex0(ijk) - gradxb(ijk)*rvvol
                  ey0(ijk) = ey0(ijk) - gradyb(ijk)*rvvol
                  ez0(ijk) = ez0(ijk) - gradzb(ijk)*rvvol
               enddo

         call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
			   dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk) - divpix(ijk)
				  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               enddo
	       write(*,*)'dopo=',dnorm



end subroutine div_cleaning

